﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace symulacja_ewolucji
{
    class Kimura
    {
        /*
         * puryny:      A G
         * piramidyny:  C T
         * 
         * tranzycje:       (alfa)
         *   puryna <-> puryna
         *   piramidyna <-> piramidyna 
         * transwersje:     (beta)
         *   puryna <-> piramidyna
         */

        /// <summary>
        /// 
        /// </summary>
        /// <param name="alfa">prawdop. zajścia tranzycji, rzędu 0.00000001</param>
        /// <param name="beta">prawdop. zajścia transwersji, rzędu 0.000000001</param>
        public Kimura(double alfa, double beta, double srednia)
        {
            this.alfa = alfa;
            this.beta = beta;
            Kimura.lambda = 1 / srednia;
        }

        private double alfa, beta;
        static double lambda;
        double[,] macierz_U = new double[4, 4];
        static Random rand = new Random();

        private void oblicz_macierz(int t, double alfa, double beta)
        {
            double s, u, r;
            s = (1 - Math.Pow(Math.E, -4 * beta * t)) / 4;
            u = (1 + Math.Pow(Math.E, -4 * beta * t) - 2*Math.Pow(Math.E, -2 * (alfa + beta) * t)) / 4;
            r = 1 - 2 * s - u;
            
            macierz_U[0, 0] = r;
            macierz_U[1, 1] = r;
            macierz_U[2, 2] = r;
            macierz_U[3, 3] = r;

            macierz_U[1, 0] = s;
            macierz_U[2, 1] = s;
            macierz_U[3, 2] = s;
            macierz_U[0, 1] = s;
            macierz_U[1, 2] = s;
            macierz_U[2, 3] = s;

            macierz_U[3, 0] = s;
            macierz_U[0, 3] = s;

            macierz_U[2, 0] = u;
            macierz_U[1, 3] = u;
            macierz_U[0, 2] = u;
            macierz_U[3, 1] = u;
        }

        int get_index(char c)
        {
            switch (c)
            {
                case 'A': return 0;
                case 'G': return 1;
                case 'C': return 2;
                case 'T': return 3;
            }
            return 0;
        }

        char get_nukleotyd(int i)
        {
            return "AGCT"[i];
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="DNA"></param>
        /// <param name="t">czas ewoluowania w latach np. 100000</param>
        /// <returns></returns>
        public string ewoluuj(string DNA, int t)
        {
            string DNA_wynikowe = "";            
            oblicz_macierz(t, alfa, beta); //rzędu 10^-8 na pokolenie
            double ilosc_podstawien = 2 * t * (alfa + beta) * DNA.Length;            
            for (int i = 0; i < DNA.Length; i++)
            {
                double p = rand.NextDouble();
                double s = 0;
                char wyn = DNA[i];
                for (int j = 0; j <= 3; j++)
                {
                    if (s < p && p < s + macierz_U[get_index(DNA[i]), j])
                    {
                        wyn = get_nukleotyd(j);
                        break;
                    }
                    s += macierz_U[get_index(DNA[i]), j];                 
                }
                DNA_wynikowe += wyn;
            }
            return DNA_wynikowe;
        }

        public static int losuj_czas_gatunku()
        {
            double x = rand.NextDouble();
           
            //double f = lambda * Math.Pow(Math.E, -lambda * x); //prawdop. wykl. w pkt. x
            //double F = 1 - lambda * Math.Pow(Math.E, -lambda * x); //dystrybuanta w pkt. u
            double czas = -Math.Log(1 - x) / lambda; // zmienne losowe o rozkl. wykł.
            //http://www.statlect.com/ucdexp1.htm

            //http://www.inf.ethz.ch/personal/anmaria/teaching/Huelsenbeck.pdf
            //double czas = -Math.Log(x) / lambda;

            return (int)czas;
        }
    }
}